"""
    GrossErrorCheck

First-stage quality control performing gross error checks on observations.
This module implements range checks, innovation checks, and basic physical
consistency validation for all observation types.

# Mathematical Foundation

## Range Check
Verify observation falls within physical bounds:
```
y_min ≤ y_obs ≤ y_max
```

## Innovation Check (3-sigma rule)
Check innovation against observation error:
```
|y_obs - H(x_b)| ≤ k * σ_obs
```
where k is typically 3-5 depending on observation type.

## Normalized Innovation
```
d_norm = (y_obs - H(x_b)) / σ_obs
```

# QC Criteria

Different observation types have different physical bounds and thresholds:

**Temperature**: 150K - 350K, innovation threshold: 10K
**Pressure**: 10hPa - 1100hPa, innovation threshold: 5hPa
**Winds**: -150 m/s - 150 m/s, innovation threshold: 15 m/s
**Humidity**: 0 - 50 g/kg, innovation threshold: 5 g/kg
**Radiance**: 100K - 400K BT, innovation threshold: 20K
"""
module GrossErrorCheck

using LinearAlgebra
using Statistics
using Dates
using ...ObservationTypes
using ...GSICoreAnalysis: AbstractStateVector

# Export types and functions
export GrossErrorConfig, GrossErrorResult
export perform_gross_error_check, check_physical_bounds, check_innovation
export get_default_bounds, get_default_thresholds

"""
    GrossErrorConfig{T<:AbstractFloat}

Configuration for gross error checking.

# Fields
- `range_bounds::Dict{ObservationVariable,Tuple{T,T}}`: Min/max physical bounds per variable
- `innovation_thresholds::Dict{ObservationVariable,T}`: Innovation thresholds per variable
- `sigma_multiplier::T`: Multiplier for sigma test (default 3.0)
- `use_normalized_innovation::Bool`: Use normalized innovation for checks
- `check_temporal_consistency::Bool`: Check time consistency
- `check_vertical_consistency::Bool`: Check vertical consistency
"""
struct GrossErrorConfig{T<:AbstractFloat}
    range_bounds::Dict{ObservationVariable,Tuple{T,T}}
    innovation_thresholds::Dict{ObservationVariable,T}
    sigma_multiplier::T
    use_normalized_innovation::Bool
    check_temporal_consistency::Bool
    check_vertical_consistency::Bool

    function GrossErrorConfig{T}() where T
        # Default physical bounds for each variable type
        bounds = Dict{ObservationVariable,Tuple{T,T}}(
            # Thermodynamic variables
            TEMPERATURE => (T(150.0), T(350.0)),           # 150K to 350K
            PRESSURE => (T(1000.0), T(110000.0)),          # 10 hPa to 1100 hPa
            SURFACE_PRESSURE => (T(50000.0), T(110000.0)), # 500 hPa to 1100 hPa
            DEW_POINT => (T(150.0), T(330.0)),            # 150K to 330K
            RELATIVE_HUMIDITY => (T(0.0), T(100.0)),       # 0% to 100%
            SPECIFIC_HUMIDITY => (T(0.0), T(0.05)),        # 0 to 50 g/kg
            MIXING_RATIO => (T(0.0), T(0.06)),            # 0 to 60 g/kg

            # Dynamic variables
            WIND_U => (T(-150.0), T(150.0)),               # ±150 m/s
            WIND_V => (T(-150.0), T(150.0)),               # ±150 m/s
            WIND_SPEED => (T(0.0), T(150.0)),              # 0 to 150 m/s
            WIND_DIRECTION => (T(0.0), T(360.0)),          # 0° to 360°
            VERTICAL_VELOCITY => (T(-10.0), T(10.0)),      # ±10 m/s

            # Radiative variables
            BRIGHTNESS_TEMPERATURE => (T(100.0), T(400.0)), # 100K to 400K
            RADIANCE => (T(0.0), T(1000.0)),               # Sensor dependent

            # Hydrologic variables
            PRECIPITATION_RATE => (T(0.0), T(500.0)),      # 0 to 500 mm/hr
            CLOUD_WATER_CONTENT => (T(0.0), T(5.0)),       # 0 to 5 g/kg
            CLOUD_FRACTION => (T(0.0), T(1.0)),            # 0 to 1

            # Specialized variables
            GPS_REFRACTIVITY => (T(0.0), T(500.0)),        # 0 to 500 N-units
            GPS_BENDING_ANGLE => (T(-0.1), T(0.1)),        # ±0.1 radians
            RADAR_REFLECTIVITY => (T(-30.0), T(80.0)),     # -30 to 80 dBZ
            RADAR_VELOCITY => (T(-100.0), T(100.0)),       # ±100 m/s
        )

        # Default innovation thresholds (absolute units)
        thresholds = Dict{ObservationVariable,T}(
            TEMPERATURE => T(10.0),              # 10 K
            PRESSURE => T(500.0),                # 5 hPa
            SURFACE_PRESSURE => T(500.0),        # 5 hPa
            DEW_POINT => T(15.0),               # 15 K
            RELATIVE_HUMIDITY => T(30.0),        # 30%
            SPECIFIC_HUMIDITY => T(0.005),       # 5 g/kg
            WIND_U => T(15.0),                  # 15 m/s
            WIND_V => T(15.0),                  # 15 m/s
            WIND_SPEED => T(15.0),              # 15 m/s
            BRIGHTNESS_TEMPERATURE => T(20.0),   # 20 K
            GPS_REFRACTIVITY => T(50.0),         # 50 N-units
            GPS_BENDING_ANGLE => T(0.01),       # 0.01 radians
            RADAR_REFLECTIVITY => T(15.0),       # 15 dBZ
            RADAR_VELOCITY => T(20.0),          # 20 m/s
        )

        new{T}(
            bounds,
            thresholds,
            T(3.0),    # sigma_multiplier
            true,      # use_normalized_innovation
            true,      # check_temporal_consistency
            true       # check_vertical_consistency
        )
    end
end

GrossErrorConfig() = GrossErrorConfig{Float64}()

"""
    GrossErrorResult{T<:AbstractFloat}

Result of gross error check for a single observation.

# Fields
- `passed::Bool`: Whether observation passed gross error check
- `failure_reason::String`: Reason for failure (if failed)
- `range_check::Bool`: Passed range check
- `innovation_check::Bool`: Passed innovation check
- `innovation::T`: Computed innovation (obs - background)
- `normalized_innovation::T`: Innovation normalized by error
- `qc_flag::Int`: QC flag (0=good, 1=suspect, 2=bad)
"""
struct GrossErrorResult{T<:AbstractFloat}
    passed::Bool
    failure_reason::String
    range_check::Bool
    innovation_check::Bool
    innovation::T
    normalized_innovation::T
    qc_flag::Int
end

"""
    perform_gross_error_check(obs::AbstractObservation{T},
                             background_value::T,
                             config::GrossErrorConfig{T}) where T

Perform comprehensive gross error check on observation.

# Arguments
- `obs`: Observation to check
- `background_value`: Background (first guess) value in observation space
- `config`: Gross error check configuration

# Returns
- `GrossErrorResult{T}`: Detailed check results

# Examples
```julia
config = GrossErrorConfig{Float64}()
result = perform_gross_error_check(obs, H_xb, config)
if result.passed
    println("Observation passed gross error check")
else
    println("Rejected: ", result.failure_reason)
end
```
"""
function perform_gross_error_check(obs::AbstractObservation{T},
                                   background_value::T,
                                   config::GrossErrorConfig{T}) where T

    # Step 1: Check physical bounds
    range_check = check_physical_bounds(obs, config)

    if !range_check
        return GrossErrorResult{T}(
            false,
            "out_of_physical_bounds",
            false,
            false,
            zero(T),
            zero(T),
            2  # Bad
        )
    end

    # Step 2: Compute innovation
    innovation = obs.observed_value - background_value

    # Step 3: Check innovation against threshold
    innovation_check_passed = check_innovation(obs, innovation, config)

    if !innovation_check_passed
        # Compute normalized innovation for diagnostic
        normalized_innov = innovation / obs.observation_error.total_error

        return GrossErrorResult{T}(
            false,
            "excessive_innovation",
            true,
            false,
            innovation,
            normalized_innov,
            2  # Bad
        )
    end

    # Compute normalized innovation
    normalized_innov = innovation / obs.observation_error.total_error

    # Check sigma test
    if config.use_normalized_innovation
        if abs(normalized_innov) > config.sigma_multiplier
            return GrossErrorResult{T}(
                false,
                "failed_sigma_test",
                true,
                false,
                innovation,
                normalized_innov,
                1  # Suspect
            )
        end
    end

    # All checks passed
    return GrossErrorResult{T}(
        true,
        "",
        true,
        true,
        innovation,
        normalized_innov,
        0  # Good
    )
end

"""
    check_physical_bounds(obs::AbstractObservation{T}, config::GrossErrorConfig{T}) where T

Check if observation value is within physical bounds for its variable type.
"""
function check_physical_bounds(obs::AbstractObservation{T}, config::GrossErrorConfig{T}) where T
    # Check if bounds are defined for this variable
    if !haskey(config.range_bounds, obs.variable)
        return true  # No bounds defined, pass by default
    end

    min_bound, max_bound = config.range_bounds[obs.variable]

    # Check basic bounds
    if obs.observed_value < min_bound || obs.observed_value > max_bound
        return false
    end

    # Additional observation-type specific checks
    if obs isa RadianceObservation
        # Check satellite radiance specific bounds
        if obs.channel <= 0
            return false
        end

        if obs.scan_angle < 0 || obs.scan_angle > 90
            return false
        end

        # Cloud fraction must be [0,1]
        if obs.cloud_fraction < 0 || obs.cloud_fraction > 1
            return false
        end
    elseif obs isa RadiosondeObservation
        # Check radiosonde specific bounds
        if obs.pressure_level <= 0
            return false
        end

        if obs.time_since_launch < 0
            return false
        end

        # Ascent rate should be reasonable
        if obs.ascent_rate < -10 || obs.ascent_rate > 20  # m/s
            return false
        end
    elseif obs isa AircraftObservation
        # Check aircraft specific bounds
        if obs.flight_level < 0
            return false
        end

        if obs.air_speed < 0 || obs.air_speed > 400  # m/s (Mach 1.2)
            return false
        end

        if obs.ground_speed < 0
            return false
        end
    elseif obs isa GPSROObservation
        # Check GPS-RO specific bounds
        if obs.impact_parameter < 0
            return false
        end
    elseif obs isa RadarObservation
        # Check radar specific bounds
        if obs.range_to_target < 0 || obs.range_to_target > 500_000  # 500 km max
            return false
        end

        if obs.azimuth_angle < 0 || obs.azimuth_angle > 360
            return false
        end

        if obs.elevation_angle < -10 || obs.elevation_angle > 90
            return false
        end
    end

    return true
end

"""
    check_innovation(obs::AbstractObservation{T}, innovation::T, config::GrossErrorConfig{T}) where T

Check if innovation is within acceptable threshold.
"""
function check_innovation(obs::AbstractObservation{T}, innovation::T, config::GrossErrorConfig{T}) where T
    # Get threshold for this variable type
    if !haskey(config.innovation_thresholds, obs.variable)
        return true  # No threshold defined, pass by default
    end

    threshold = config.innovation_thresholds[obs.variable]

    # Platform-specific threshold adjustments
    threshold = adjust_threshold_for_platform(obs, threshold, config)

    # Check absolute innovation
    if abs(innovation) > threshold
        return false
    end

    return true
end

"""
    adjust_threshold_for_platform(obs::AbstractObservation{T}, threshold::T, config::GrossErrorConfig{T}) where T

Adjust innovation threshold based on platform type and conditions.
"""
function adjust_threshold_for_platform(obs::AbstractObservation{T}, threshold::T, config::GrossErrorConfig{T}) where T
    adjusted = threshold

    # Satellite observations may need relaxed thresholds
    if obs isa RadianceObservation
        # Relax threshold for high scan angles (limb observations)
        if obs.scan_angle > 45.0
            adjusted *= T(1.5)
        end

        # Relax threshold for cloudy scenes
        if obs.cloud_fraction > 0.3
            adjusted *= T(2.0)
        end
    end

    # Radiosonde observations
    if obs isa RadiosondeObservation
        # Tighter threshold for mandatory levels
        if obs.level_type == :mandatory
            adjusted *= T(0.8)
        end

        # Relax threshold for high altitude (lower pressure)
        if obs.pressure_level < 10000  # < 100 hPa
            adjusted *= T(1.3)
        end
    end

    # Aircraft observations
    if obs isa AircraftObservation
        # Different thresholds for different flight phases
        if obs.flight_phase == :climb || obs.flight_phase == :descent
            adjusted *= T(1.2)  # More variability during ascent/descent
        end
    end

    # GPS-RO observations
    if obs isa GPSROObservation
        # Relax threshold near surface (high refractivity)
        if obs.variable == GPS_REFRACTIVITY && obs.location.altitude < 5000
            adjusted *= T(1.5)
        end
    end

    return adjusted
end

"""
    get_default_bounds(var::ObservationVariable) -> Tuple{Float64, Float64}

Get default physical bounds for an observation variable.
"""
function get_default_bounds(var::ObservationVariable)
    config = GrossErrorConfig{Float64}()
    if haskey(config.range_bounds, var)
        return config.range_bounds[var]
    else
        return (-Inf, Inf)
    end
end

"""
    get_default_thresholds(var::ObservationVariable) -> Float64

Get default innovation threshold for an observation variable.
"""
function get_default_thresholds(var::ObservationVariable)
    config = GrossErrorConfig{Float64}()
    if haskey(config.innovation_thresholds, var)
        return config.innovation_thresholds[var]
    else
        return Inf
    end
end

"""
    check_temporal_consistency(obs::AbstractObservation{T}, time_window::T) where T

Check if observation time is within acceptable time window.
"""
function check_temporal_consistency(obs::AbstractObservation{T},
                                   analysis_time::DateTime,
                                   time_window::T) where T
    time_diff = abs(Dates.value(obs.observation_time - analysis_time)) / 1000.0  # seconds
    return time_diff <= time_window
end

"""
    check_vertical_consistency(obs::AbstractObservation{T}) where T

Check vertical consistency for upper-air observations.
"""
function check_vertical_consistency(obs::AbstractObservation{T}) where T
    # For radiosonde, check pressure-altitude consistency
    if obs isa RadiosondeObservation
        # Approximate standard atmosphere
        # p = p0 * exp(-z/H) where H ≈ 7400 m
        expected_pressure = 101325.0 * exp(-obs.location.altitude / 7400.0)

        # Allow 20% deviation
        if abs(obs.pressure_level - expected_pressure) / expected_pressure > 0.2
            return false
        end
    end

    return true
end

end # module GrossErrorCheck
